目次

# conda install -c conda-forge osmnx

# https://osmnx.readthedocs.io/en/stable/
# https://qiita.com/moshi/items/e383491664d028cd2166
# 
# https://zenn.dev/homing/articles/f9a314841c737d
# https://qiita.com/naozo-se/items/bc168a2282907a6b29bf


# https://stackoverflow.com/questions/62285134/how-to-fill-water-bodies-with-osmnx-in-python
import osmnx as ox
import networkx as nx
import folium
from geopy.distance import distance

# 与えられた緯度経度

latitude = 34.69155489576708  # 例としてケーニッヒスベルクの緯度
longitude = 135.48959938567188  # 例としてケーニッヒスベルクの経度

# 2kmの範囲を定義
distance_km = 2

# 地図データを取得
G = ox.graph_from_point((latitude, longitude),
                        dist=distance_km*1000)
   #network_type='drive')

# グラフを取得
nodes, edges = ox.graph_to_gdfs(G)

# 橋の情報を抽出('bridge'タグを持つエッジのみを抽出)
bridge_edges = edges[edges['bridge'] == 'yes']
#bridge_edges = edges

# グラフのノードとエッジを抽出
bridge_graph = nx.Graph()

# エッジ(橋)をグラフに追加
for idx, data in bridge_edges.iterrows():
    u, v, key = idx
    bridge_graph.add_edge(u, v)

# オイラー閉路の判定
eulerian = nx.is_eulerian(bridge_graph)

# Foliumで可視化
m = folium.Map(location=[latitude, longitude], zoom_start=14)

# グラフのエッジ(橋)を地図に追加
for u, v in bridge_graph.edges():
    point_u = (G.nodes[u]['y'], G.nodes[u]['x'])
    point_v = (G.nodes[v]['y'], G.nodes[v]['x'])
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(m)

# 中心地点にマーカーを追加
folium.Marker(location=[latitude, longitude], popup="中心地点").add_to(m)

# オイラー閉路の有無を表示
if eulerian:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路が存在します。", icon=folium.Icon(color='green')).add_to(m)
else:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路は存在しません。", icon=folium.Icon(color='red')).add_to(m)

m
Make this Notebook Trusted to load map: File -> Trust Notebook
import osmnx as ox
import networkx as nx
import folium

# 与えられた緯度経度
latitude = 54.7088  # 例としてケーニッヒスベルクの緯度
longitude = 20.5108  # 例としてケーニッヒスベルクの経度

# 2kmの範囲を定義
distance_km = 2

# 地図データを取得
G = ox.graph_from_point((latitude, longitude), dist=distance_km*1000, network_type='all')

# 無向グラフに変換
G = G.to_undirected()

# グラフの連結成分を取得
connected_components = list(nx.connected_components(G))

# 橋の情報を抽出('bridge'タグを持つエッジのみを抽出)
bridge_edges = [edge for edge in G.edges(data=True) if 'bridge' in edge[2] and edge[2]['bridge'] == 'yes']

# 各連結成分を一つのノード(陸地)と見なす
land_nodes = {i: component for i, component in enumerate(connected_components)}

# 橋のノードを含む陸地ノードのIDを見つける関数
def find_land_node(node):
    for land_id, nodes in land_nodes.items():
        if node in nodes:
            return land_id
    return None

# 橋で陸地ノードを接続するための新しいグラフを作成
final_graph = nx.Graph()

# 陸地ノードを追加
for land_id, nodes in land_nodes.items():
    final_graph.add_node(land_id, nodes=nodes)

# 橋エッジを追加(橋が異なる陸地を接続している場合)
for u, v, data in bridge_edges:
    land_u = find_land_node(u)
    land_v = find_land_node(v)
    if land_u is not None and land_v is not None and land_u != land_v:
        final_graph.add_edge(land_u, land_v, **data)

# オイラー閉路の判定
eulerian = nx.is_eulerian(final_graph)

# Foliumで可視化
m = folium.Map(location=[latitude, longitude], zoom_start=14)

# 橋のエッジを地図に追加
for u, v, data in final_graph.edges(data=True):
    land_u_nodes = land_nodes[u]
    land_v_nodes = land_nodes[v]
    point_u = next((G.nodes[node]['y'], G.nodes[node]['x']) for node in land_u_nodes)
    point_v = next((G.nodes[node]['y'], G.nodes[node]['x']) for node in land_v_nodes)
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(m)

# 中心地点にマーカーを追加
folium.Marker(location=[latitude, longitude], popup="中心地点").add_to(m)

# オイラー閉路の有無を表示
if eulerian:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路が存在します。", icon=folium.Icon(color='green')).add_to(m)
else:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路は存在しません。", icon=folium.Icon(color='red')).add_to(m)

m
Make this Notebook Trusted to load map: File -> Trust Notebook
print(land_graph)
print(final_graph)
MultiGraph with 5819 nodes and 8690 edges
Graph with 1 nodes and 0 edges
import folium
import osmnx as ox

import geopandas as gpd

# 対象地域の道路情報取得 (愛知県名古屋市中村区)
#query = "Nakamuraku,Nagoya,Aichi,Japan"
#G = ox.graph_from_place(query, network_type="drive")

(lat,lon)=(34.69199639822432, 135.49084584150586)

# Bounding box
bN,bS,bE,bW=lon+1,lon-1,lon-1,lon+1

#graph=ox.graph.graph_from_point((lat, lon),
#                                dist=1000,
#                                network_type='drive')
graph=ox.geometries.geometries_from_bbox(bN, bS, bE, bW,
                                tags={"natural": "water"})

# 道路グラフネットワーク可視化
#fmap=ox.plot_graph_folium(graph)
fmap=ox.plot_graph_folium(graph)
folium.LayerControl().add_to(fmap)
fmap


#fmap.save(outfile="road_network.html")
#opts = {"node_size": 5, "bgcolor": "white", "node_color": "blue", "edge_color": "blue"}
#ox.plot_graph(G, show=False, save=True, filepath="road_network.png", **opts)
/var/folders/mn/nmxx09wj21j7f81mwghvt1sw0000gn/T/ipykernel_40987/1645265506.py:18: FutureWarning: The `geometries` module and `geometries_from_X` functions have been renamed the `features` module and `features_from_X` functions. Use these instead. The `geometries` module and function names are deprecated and will be removed in the v2.0.0 release. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  graph=ox.geometries.geometries_from_bbox(bN, bS, bE, bW,
/opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/geometries.py:48: FutureWarning: The `north`, `south`, `east`, and `west` parameters are deprecated and will be removed in the v2.0.0 release. Use the `bbox` parameter instead. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  return features.features_from_bbox(north, south, east, west, tags=tags)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/geopandas/array.py:950, in GeometryArray.estimate_utm_crs(self, datum_name)
    949 try:
--> 950     return CRS.from_epsg(utm_crs_list[0].code)
    951 except IndexError:

IndexError: list index out of range

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
Cell In[17], line 18
     13 bN,bS,bE,bW=lon+1,lon-1,lon-1,lon+1
     15 #graph=ox.graph.graph_from_point((lat, lon),
     16 #                                dist=1000,
     17 #                                network_type='drive')
---> 18 graph=ox.geometries.geometries_from_bbox(bN, bS, bE, bW,
     19                                 tags={"natural": "water"})
     21 # 道路グラフネットワーク可視化
     22 #fmap=ox.plot_graph_folium(graph)
     23 fmap=ox.plot_graph_folium(graph)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/geometries.py:48, in geometries_from_bbox(north, south, east, west, tags)
     22 """
     23 Do not use: deprecated.
     24 
   (...)
     45 gdf : geopandas.GeoDataFrame
     46 """
     47 warn(DEP_MSG, FutureWarning, stacklevel=2)
---> 48 return features.features_from_bbox(north, south, east, west, tags=tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:134, in features_from_bbox(north, south, east, west, bbox, tags)
    131 polygon = utils_geo.bbox_to_poly(bbox=bbox)
    133 # create GeoDataFrame of features within this polygon
--> 134 return features_from_polygon(polygon, tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:349, in features_from_polygon(polygon, tags)
    346 response_jsons = _overpass._download_overpass_features(polygon, tags)
    348 # create GeoDataFrame from the downloaded data
--> 349 return _create_gdf(response_jsons, polygon, tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:440, in _create_gdf(response_jsons, polygon, tags)
    437 relation_types = {"boundary", "multipolygon"}
    439 # extract geometries from the downloaded osm data
--> 440 for response_json in response_jsons:
    441     response_count += 1
    443     # Parses the JSON of OSM nodes, ways and (multipolygon) relations
    444     # to dictionaries of coordinates, Shapely Points, LineStrings,
    445     # Polygons and MultiPolygons

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/_overpass.py:389, in _download_overpass_features(polygon, tags)
    373 """
    374 Retrieve OSM features within boundary from the Overpass API.
    375 
   (...)
    386     a generator of JSON responses from the Overpass server
    387 """
    388 # subdivide query polygon to get list of sub-divided polygon coord strings
--> 389 polygon_coord_strs = _make_overpass_polygon_coord_strs(polygon)
    390 utils.log(f"Requesting data from API in {len(polygon_coord_strs)} request(s)")
    392 # pass exterior coordinates of each polygon in list to API, one at a time

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/_overpass.py:253, in _make_overpass_polygon_coord_strs(polygon)
    234 """
    235 Subdivide query polygon and return list of coordinate strings.
    236 
   (...)
    249     list of strings of exterior coordinates of polygon(s)
    250 """
    251 # first subdivide the polygon if its area exceeds max size
    252 # this results in a multipolygon of 1+ constituent polygons
--> 253 poly_proj, crs_proj = projection.project_geometry(polygon)
    254 multi_poly_proj = utils_geo._consolidate_subdivide_geometry(poly_proj)
    255 multi_poly, _ = projection.project_geometry(multi_poly_proj, crs=crs_proj, to_latlong=True)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/projection.py:60, in project_geometry(geometry, crs, to_crs, to_latlong)
     57     crs = settings.default_crs
     59 gdf = gpd.GeoDataFrame(geometry=[geometry], crs=crs)
---> 60 gdf_proj = project_gdf(gdf, to_crs=to_crs, to_latlong=to_latlong)
     61 geometry_proj = gdf_proj["geometry"].iloc[0]
     62 return geometry_proj, gdf_proj.crs

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/projection.py:99, in project_gdf(gdf, to_crs, to_latlong)
     97 # else if to_crs is None, project gdf to an appropriate UTM zone
     98 elif to_crs is None:
---> 99     to_crs = gdf.estimate_utm_crs()
    101 # project the gdf
    102 gdf_proj = gdf.to_crs(to_crs)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/geopandas/geodataframe.py:1451, in GeoDataFrame.estimate_utm_crs(self, datum_name)
   1415 def estimate_utm_crs(self, datum_name="WGS 84"):
   1416     """Returns the estimated UTM CRS based on the bounds of the dataset.
   1417 
   1418     .. versionadded:: 0.9
   (...)
   1449     - Prime Meridian: Greenwich
   1450     """
-> 1451     return self.geometry.estimate_utm_crs(datum_name=datum_name)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/geopandas/geoseries.py:1214, in GeoSeries.estimate_utm_crs(self, datum_name)
   1178 def estimate_utm_crs(self, datum_name: str = "WGS 84") -> CRS:
   1179     """Returns the estimated UTM CRS based on the bounds of the dataset.
   1180 
   1181     .. versionadded:: 0.9
   (...)
   1212     - Prime Meridian: Greenwich
   1213     """
-> 1214     return self.values.estimate_utm_crs(datum_name)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/geopandas/array.py:952, in GeometryArray.estimate_utm_crs(self, datum_name)
    950     return CRS.from_epsg(utm_crs_list[0].code)
    951 except IndexError:
--> 952     raise RuntimeError("Unable to determine UTM CRS")

RuntimeError: Unable to determine UTM CRS
import osmnx as ox
ox.settings.log_console = True

# add more items here to get all the landforms you want
places = ['Manhattan, NY, USA', 'Brooklyn, NY, USA', 'Queens, NY, USA', 'Bronx, NY, USA']
land = ox.geocode_to_gdf(places)

# get the water bodies
left, bottom, right, top = land.total_bounds
bbox = top, bottom, right, left
poly = ox.utils_geo.bbox_to_poly(*bbox)
water = ox.features_from_polygon(poly, tags={'natural': 'water'})

# constrain the plotting window as desired
c = land.unary_union.centroid
bbox = ox.utils_geo.bbox_from_point((c.y, c.x), dist=12000)

water_color = 'blue'
land_color = '#aaaaaa'
fig, ax = ox.plot_footprints(water, bbox=bbox,
                             color=water_color, bgcolor=water_color,
                             show=False, close=False)
ax = land.plot(ax=ax, zorder=0, fc=land_color)
/var/folders/mn/nmxx09wj21j7f81mwghvt1sw0000gn/T/ipykernel_40987/273932149.py:11: FutureWarning: The `north`, `south`, `east`, and `west` parameters are deprecated and will be removed in the v2.0.0 release. Use the `bbox` parameter instead. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  poly = ox.utils_geo.bbox_to_poly(*bbox)
_images/e83ae0b43aa451f42a7e88bffda1d8bb9cdecb503c883c37fea9cd9735049e04.png
import osmnx as ox
import folium

# 場所を設定(例:東京都渋谷区)
place_name = "Shibuya, Tokyo, Japan"

(lat,lon)=(34.69199639822432, 135.49084584150586)

# 陸地と道路のグラフを取得
#G = ox.graph_from_place(place_name, network_type='drive')
G = ox.graph_from_point((lat,lon), network_type='drive')

nodes, edges = ox.graph_to_gdfs(G)

# 中心点を取得
lat = nodes['y'].mean()
lon = nodes['x'].mean()

# foliumマップを作成
m = folium.Map(location=[lat, lon], zoom_start=14)

# 陸地を描画
folium.GeoJson(data=edges).add_to(m)

# マップを保存または表示
m.save('map.html')

# マップを表示(Jupyter Notebookの場合)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
import osmnx as ox
import folium
import geopandas as gpd

# 場所を設定(例:東京都渋谷区)
place_name = "Shibuya, Tokyo, Japan"

# 陸地のポリゴンを取得
land_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': 'land'})

# 水域のポリゴンを取得
water_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': ['water', 'bay', 'strait']})

# 道路のグラフを取得
G = ox.graph_from_place(place_name, network_type='drive')
nodes, edges = ox.graph_to_gdfs(G)

# 中心点を取得
lat = nodes['y'].mean()
lon = nodes['x'].mean()

# foliumマップを作成
m = folium.Map(location=[lat, lon], zoom_start=14)

# 陸地を描画
for _, polygon in land_polygons.iterrows():
    folium.GeoJson(polygon['geometry'], style_function=lambda x: {'fillColor': 'green', 'color': 'green'}).add_to(m)

# 水域を描画
for _, polygon in water_polygons.iterrows():
    folium.GeoJson(polygon['geometry'], style_function=lambda x: {'fillColor': 'blue', 'color': 'blue'}).add_to(m)

# 道路を描画
folium.GeoJson(data=edges).add_to(m)

# マップを保存または表示
m.save('map.html')

# マップを表示(Jupyter Notebookの場合)
m
/var/folders/mn/nmxx09wj21j7f81mwghvt1sw0000gn/T/ipykernel_40987/388463500.py:9: FutureWarning: The `geometries` module and `geometries_from_X` functions have been renamed the `features` module and `features_from_X` functions. Use these instead. The `geometries` module and function names are deprecated and will be removed in the v2.0.0 release. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  land_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': 'land'})
---------------------------------------------------------------------------
InsufficientResponseError                 Traceback (most recent call last)
Cell In[22], line 9
      6 place_name = "Shibuya, Tokyo, Japan"
      8 # 陸地のポリゴンを取得
----> 9 land_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': 'land'})
     11 # 水域のポリゴンを取得
     12 water_polygons = ox.geometries.geometries_from_place(place_name, tags={'natural': ['water', 'bay', 'strait']})

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/geometries.py:128, in geometries_from_place(query, tags, which_result, buffer_dist)
    104 """
    105 Do not use: deprecated.
    106 
   (...)
    125 gdf : geopandas.GeoDataFrame
    126 """
    127 warn(DEP_MSG, FutureWarning, stacklevel=2)
--> 128 return features.features_from_place(query, tags, which_result, buffer_dist)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:298, in features_from_place(query, tags, which_result, buffer_dist)
    295 utils.log("Constructed place geometry polygon(s) to query API")
    297 # create GeoDataFrame using this polygon(s) geometry
--> 298 return features_from_polygon(polygon, tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:349, in features_from_polygon(polygon, tags)
    346 response_jsons = _overpass._download_overpass_features(polygon, tags)
    348 # create GeoDataFrame from the downloaded data
--> 349 return _create_gdf(response_jsons, polygon, tags)

File /opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/features.py:489, in _create_gdf(response_jsons, polygon, tags)
    487 if len(geometries) == 0:  # pragma: no cover
    488     msg = "No data elements in server response. Check log and query location/tags."
--> 489     raise InsufficientResponseError(msg)
    491 # remove untagged elements from the final dict of geometries
    492 utils.log(f"{len(geometries)} geometries created in the dict")

InsufficientResponseError: No data elements in server response. Check log and query location/tags.
# 試しコード
import osmnx as ox
import networkx as nx
import folium

# 地図データを取得
# c.f. https://osmnx.readthedocs.io/en/stable/user-reference.html#module-osmnx.graph
graph=ox.graph.graph_from_bbox(
    # bbox : (north, south, east, west)
    bbox=(35.05,34.90,
          135.75,135.80), 
    # "all", "all_public", "bike", "drive", "drive_service", "walk"
    network_type='bike',
    simplify=True, # True
    retain_all=False,
    truncate_by_edge=False,
    custom_filter=None)
# 無向グラフにしておく
graph=graph.to_undirected()
# グラフを取得
nodes, edges = ox.graph_to_gdfs(graph)
# グラフのノードとエッジを追加
nxgraph = nx.Graph()
for idx, data in edges.iterrows():
    u, v, key = idx
    nxgraph.add_edge(u, v)
    
# Foliumで可視化
lat, lon = [(35.05+34.90)/2,
                       (135.80+135.75)/2]
fmap = folium.Map(location=[latitude, longitude], zoom_start=13)
# グラフを地図に追加
for u, v in nxgraph.edges():
    point_u = (graph.nodes[u]['y'],graph.nodes[u]['x'])
    point_v = (graph.nodes[v]['y'],graph.nodes[v]['x'])
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(fmap)
fmap
Make this Notebook Trusted to load map: File -> Trust Notebook
# 橋に絞る
bridge_edges = edges[edges['bridge'] == 'yes']
# グラフのノードとエッジを追加
bridge_nxgraph = nx.Graph()
for idx, data in bridge_edges.iterrows():
    u, v, key = idx
    bridge_nxgraph.add_edge(u, v)
    
# Foliumで可視化
lat, lon = [(35.05+34.90)/2,
                       (135.80+135.75)/2]
#attr=('&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> '
#    'contributors, &copy; <a href="https://cartodb.com/attributions">CartoDB</a>')
#tiles = "https://{s}.basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}.png"
#fmap = folium.Map(location=[lat, lon], tiles=tiles, attr=attr, zoom_start=13)
fmap = folium.Map(location=[latitude, longitude], zoom_start=13)
# グラフを地図に追加
for u, v in nxgraph.edges():
    point_u = (graph.nodes[u]['y'],graph.nodes[u]['x'])
    point_v = (graph.nodes[v]['y'],graph.nodes[v]['x'])
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(fmap)
# グラフを地図に追加
for u, v in bridge_nxgraph.edges():
    point_u = (graph.nodes[u]['y'],graph.nodes[u]['x'])
    point_v = (graph.nodes[v]['y'],graph.nodes[v]['x'])
    folium.PolyLine(locations=[point_u, point_v], color='red').add_to(fmap)
fmap
Make this Notebook Trusted to load map: File -> Trust Notebook
# 橋の情報を抽出('bridge'タグを持つエッジのみを抽出)
bridge_edges2 = [(u, v) for u, v, data in graph.edges(data=True) \
                if 'bridge' in data and data['bridge'] == 'yes']

# 橋のノードを収集
bridge_nodes = set()
for u, v in bridge_edges2:
    bridge_nodes.add(u)
    bridge_nodes.add(v)

# 橋を除いたグラフを作成
land_graph = graph.copy()
land_graph.remove_edges_from(bridge_edges2)


# 陸地の連結成分を取得
land_components = list(nx.connected_components(land_graph))

# 各連結成分を一つのノード(陸地)と見なす
land_nodes = {i: component for i, component in enumerate(land_components)}

# 橋で陸地ノードを接続するための新しいグラフを作成
final_graph = nx.Graph()
# 陸地ノードを追加
for land_id, nodes in land_nodes.items():
    final_graph.add_node(land_id, nodes=nodes)
# 橋エッジを追加(橋が異なる陸地を接続している場合)
for u, v in bridge_edges2:
    land_u = next((land_id for land_id, nodes in land_nodes.items() if u in nodes), None)
    land_v = next((land_id for land_id, nodes in land_nodes.items() if v in nodes), None)
    if land_u is not None and land_v is not None and land_u != land_v:
        final_graph.add_edge(land_u, land_v)
        #print('added')

# Foliumで可視化
m = folium.Map(location=[latitude, longitude], zoom_start=12)
# 橋のエッジを地図に追加
for u, v in final_graph.edges():
    land_u_nodes = land_nodes[u]
    land_v_nodes = land_nodes[v]
    point_u = next((graph.nodes[node]['y'],graph.nodes[node]['x']) for node in land_u_nodes)
    point_v = next((graph.nodes[node]['y'],graph.nodes[node]['x']) for node in land_v_nodes)
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
import networkx as nx
import matplotlib.pyplot as plt

# final_graphの可視化
plt.figure(figsize=(10, 8))
pos = nx.spring_layout(final_graph)  # ノードの位置を計算

# ノードを描画
nx.draw_networkx_nodes(final_graph, pos, node_size=300, node_color='lightblue')

# エッジを描画
nx.draw_networkx_edges(final_graph, pos, width=2, edge_color='blue')

# ノードのラベルを描画
nx.draw_networkx_labels(final_graph, pos, font_size=12, font_family='sans-serif')

# グラフを表示
plt.title("Final Graph Visualization")
plt.axis('off')
plt.show()
_images/e602ca3b7740976eb4ae47afefe516f29376c84dc3657fd3ffa523aca229d631.png
import osmnx as ox

# Get a river network and plot it with all edge intersections.
point = (35.06+34.99)/2,(135.79+135.75)/2 # lat, long
G = ox.graph_from_point(point, 
                        dist=6000,
                        custom_filter='["waterway"~"river"]')
fig, ax = ox.plot_graph(G, node_color='r')
fig, ax = ox.plot_graph(G, node_color='r')
_images/d5643776b3d4427b46f9aa3a70caeea29848be9762ec6f7db13f70cd453c1bac.png
import osmnx as ox

# Get a river network and plot it with all edge intersections.
point = (35.06+34.99)/2,(135.79+135.75)/2 # lat, long
G = ox.graph_from_point(point, 
                        dist=10000,
                        custom_filter='["waterway"~"river"]')
fig, ax = ox.plot_graph(G, node_color='r')
_images/4ca70e5142721b24c5ba8217dfed324c2c93b728904d1bac32161b7011ac0a23.png
G2 = ox.graph_from_point(point, 
                        dist=200000,
                        custom_filter='["bridge"]')
print(G2)
fig, ax = ox.plot_graph(G2, node_color='r')
/opt/anaconda3/envs/jupyter-book-py311/lib/python3.11/site-packages/osmnx/_overpass.py:254: UserWarning: This area is 64 times your configured Overpass max query area size. It will automatically be divided up into multiple sub-queries accordingly. This may take a long time.
  multi_poly_proj = utils_geo._consolidate_subdivide_geometry(poly_proj)
import osmnx as ox
ox.settings.log_console = True

# add more items here to get all the landforms you want
places = ['左京区, 京都市, 京都府, Japan']

land = ox.geocode_to_gdf(places)

# get the water bodies
left, bottom, right, top = land.total_bounds
bbox = top, bottom, right, left
poly = ox.utils_geo.bbox_to_poly(*bbox)
water = ox.features_from_polygon(poly, tags={'natural': 'water'})

# constrain the plotting window as desired
c = land.unary_union.centroid
bbox = ox.utils_geo.bbox_from_point((c.y, c.x), dist=12000)

water_color = 'blue'
land_color = '#aaaaaa'
fig, ax = ox.plot_footprints(water, bbox=bbox,
                             color=water_color, bgcolor=water_color,
                             show=False, close=False)
ax = land.plot(ax=ax, zorder=0, fc=land_color)
/var/folders/mn/nmxx09wj21j7f81mwghvt1sw0000gn/T/ipykernel_1063/3265604966.py:12: FutureWarning: The `north`, `south`, `east`, and `west` parameters are deprecated and will be removed in the v2.0.0 release. Use the `bbox` parameter instead. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  poly = ox.utils_geo.bbox_to_poly(*bbox)
_images/39e6be15c6a10bbc7cc4b72ba67e21c011bcd82e7e3d851878e2a34bde19bfc1.png
# 橋の情報を抽出('bridge'タグを持つエッジのみを抽出)
bridge_edges = [(u, v) for u, v, data in graph.edges(data=True) \
                if 'bridge' in data and data['bridge'] == 'yes']

# 橋のノードを収集
bridge_nodes = set()
for u, v in bridge_edges:
    bridge_nodes.add(u)
    bridge_nodes.add(v)

# 橋を除いたグラフを作成
land_graph = graph.copy()
land_graph.remove_edges_from(bridge_edges)
# 陸地の連結成分を取得
land_components = list(nx.connected_components(land_graph))
# 各連結成分を一つのノード(陸地)と見なす
land_nodes = {i: component for i, component in enumerate(land_components)}
# 橋で陸地ノードを接続するための新しいグラフを作成
final_graph = nx.Graph()
# 陸地ノードを追加
for land_id, nodes in land_nodes.items():
    final_graph.add_node(land_id, nodes=nodes)
# 橋エッジを追加(橋が異なる陸地を接続している場合)
for u, v in bridge_edges:
    land_u = next((land_id for land_id, nodes in land_nodes.items() if u in nodes), None)
    land_v = next((land_id for land_id, nodes in land_nodes.items() if v in nodes), None)
    if land_u is not None and land_v is not None and land_u != land_v:
        final_graph.add_edge(land_u, land_v)

# Foliumで可視化
m = folium.Map(location=[latitude, longitude], zoom_start=12)
# 橋のエッジを地図に追加
for u, v in final_graph.edges():
    land_u_nodes = land_nodes[u]
    land_v_nodes = land_nodes[v]
    point_u = next((graph.nodes[node]['y'],graph.nodes[node]['x']) for node in land_u_nodes)
    point_v = next((graph.nodes[node]['y'],graph.nodes[node]['x']) for node in land_v_nodes)
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
# 橋の情報を抽出('bridge'タグを持つエッジのみを抽出)
#bridge_edges = edges[edges['bridge'] == 'yes']
bridge_edges = edges

# グラフのノードとエッジを抽出
bridge_graph = nx.Graph()

# エッジ(橋)をグラフに追加
for idx, data in bridge_edges.iterrows():
    u, v, key = idx
    bridge_graph.add_edge(u, v)

# オイラー閉路の判定
eulerian = nx.is_eulerian(bridge_graph)
import osmnx as ox
import networkx as nx
import folium

# 与えられた緯度経度
latitude = 54.7088  # 例としてケーニッヒスベルクの緯度
longitude = 20.5108  # 例としてケーニッヒスベルクの経度

# 2kmの範囲を定義
distance_km = 2

# 地図データを取得
# c.f. https://osmnx.readthedocs.io/en/stable/user-reference.html#module-osmnx.graph
G = ox.graph_from_point((latitude, longitude), dist=distance_km*1000, network_type='all')

# 無向グラフに変換
G = G.to_undirected()

# 橋の情報を抽出('bridge'タグを持つエッジのみを抽出)
bridge_edges = [(u, v) for u, v, data in G.edges(data=True) if 'bridge' in data and data['bridge'] == 'yes']

# 橋のノードを収集
bridge_nodes = set()
for u, v in bridge_edges:
    bridge_nodes.add(u)
    bridge_nodes.add(v)

# 橋を除いたグラフを作成
land_graph = G.copy()
land_graph.remove_edges_from(bridge_edges)

# 陸地の連結成分を取得
land_components = list(nx.connected_components(land_graph))

# 各連結成分を一つのノード(陸地)と見なす
land_nodes = {i: component for i, component in enumerate(land_components)}

# 橋で陸地ノードを接続するための新しいグラフを作成
final_graph = nx.Graph()

# 陸地ノードを追加
for land_id, nodes in land_nodes.items():
    final_graph.add_node(land_id, nodes=nodes)

# 橋エッジを追加(橋が異なる陸地を接続している場合)
for u, v in bridge_edges:
    land_u = next((land_id for land_id, nodes in land_nodes.items() if u in nodes), None)
    land_v = next((land_id for land_id, nodes in land_nodes.items() if v in nodes), None)
    if land_u is not None and land_v is not None and land_u != land_v:
        final_graph.add_edge(land_u, land_v)

# オイラー閉路の判定
eulerian = nx.is_eulerian(final_graph)

# Foliumで可視化
m = folium.Map(location=[latitude, longitude], zoom_start=14)

# 橋のエッジを地図に追加
for u, v in final_graph.edges():
    land_u_nodes = land_nodes[u]
    land_v_nodes = land_nodes[v]
    point_u = next((G.nodes[node]['y'], G.nodes[node]['x']) for node in land_u_nodes)
    point_v = next((G.nodes[node]['y'], G.nodes[node]['x']) for node in land_v_nodes)
    folium.PolyLine(locations=[point_u, point_v], color='blue').add_to(m)

# 中心地点にマーカーを追加
folium.Marker(location=[latitude, longitude], popup="中心地点").add_to(m)

# オイラー閉路の有無を表示
if eulerian:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路が存在します。", icon=folium.Icon(color='green')).add_to(m)
else:
    folium.Marker(location=[latitude, longitude], popup="オイラー閉路は存在しません。", icon=folium.Icon(color='red')).add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
import networkx as nx
import matplotlib.pyplot as plt

# final_graphの可視化
plt.figure(figsize=(10, 8))
pos = nx.spring_layout(final_graph)  # ノードの位置を計算

# ノードを描画
nx.draw_networkx_nodes(final_graph, pos, node_size=300, node_color='lightblue')

# エッジを描画
nx.draw_networkx_edges(final_graph, pos, width=2, edge_color='blue')

# ノードのラベルを描画
nx.draw_networkx_labels(final_graph, pos, font_size=12, font_family='sans-serif')

# グラフを表示
plt.title("Final Graph Visualization")
plt.axis('off')
plt.show()
_images/405543c51cc32693a192b3fa85b01bb5743a00b1a2d6356e56f8e9c2f1be3dd5.png